Cross species correlation

Notebook for cell type identity correlations

General Setup

Setup chunk

Setup reticulate

knitr::opts_chunk$set(fig.width = 8)
knitr::opts_knit$set(root.dir = normalizePath(".."))
knitr::opts_knit$get("root.dir")

[1] "/nas/groups/treutlein/USERS/tomasgomes/projects/pallium_evo"

Load libraries

library(reticulate)
knitr::knit_engines$set(python = reticulate::eng_python)
py_available(initialize = FALSE)

[1] FALSE

use_python(Sys.which("python"))
py_config()

python:         /home/tpires/bin/miniconda3/bin/python
libpython:      /home/tpires/bin/miniconda3/lib/libpython3.8.so
pythonhome:     /home/tpires/bin/miniconda3:/home/tpires/bin/miniconda3
version:        3.8.3 (default, May 19 2020, 18:47:26)  [GCC 7.3.0]
numpy:          /home/tpires/bin/miniconda3/lib/python3.8/site-packages/numpy
numpy_version:  1.20.2

NOTE: Python version was forced by RETICULATE_PYTHON

Load and reformat data

Load Seurat data

library(igraph)

Attaching package: ‘igraph’

The following objects are masked from ‘package:purrr’:

    compose, simplify

The following object is masked from ‘package:tidyr’:

    crossing

The following object is masked from ‘package:tibble’:

    as_data_frame

The following objects are masked from ‘package:dplyr’:

    as_data_frame, groups, union

The following objects are masked from ‘package:stats’:

    decompose, spectrum

The following object is masked from ‘package:base’:

    union

Load axolotl data

axolotl_gaba = readRDS(file = "data/expression/axolotl/pallium_neuronal_GABA_res07_harmony.RDS")
axolotl_glut = readRDS(file = "data/expression/axolotl/pallium_neuronal_Glut_res05_harmony.RDS")
axolotl_neurons = readRDS(file = "data/expression/axolotl/pallium_neuronal_res07_harmony.RDS")
axolotl_all = readRDS(file = "data/expression/axolotl/palliumres07_harmony.RDS")

Merge axolotl annotations together into a single object

```r
lizard_all = readRDS(\data/expression/lizard_all_v3.RDS\)
turtle_all = readRDS(\data/expression/turtle_all_v3.RDS\)
drerio_brain = readRDS(\data/expression/drerio_brain_v3.RDS\)
finches_brain = readRDS(\data/expression/HVC_RA_X.RDS\)
mouse_brain = readRDS(\data/expression/l5_all_seurat.RDS\)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Organise metadata into a list


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuYXhvbG90bF9nYWJhID0gcmVhZFJEUyhmaWxlID0gXFxkYXRhL2V4cHJlc3Npb24vYXhvbG90bC9wYWxsaXVtX25ldXJvbmFsX0dBQkFfcmVzMDdfaGFybW9ueS5SRFNcXClcbmF4b2xvdGxfZ2x1dCA9IHJlYWRSRFMoZmlsZSA9IFxcZGF0YS9leHByZXNzaW9uL2F4b2xvdGwvcGFsbGl1bV9uZXVyb25hbF9HbHV0X3JlczA1X2hhcm1vbnkuUkRTXFwpXG5heG9sb3RsX25ldXJvbnMgPSByZWFkUkRTKGZpbGUgPSBcXGRhdGEvZXhwcmVzc2lvbi9heG9sb3RsL3BhbGxpdW1fbmV1cm9uYWxfcmVzMDdfaGFybW9ueS5SRFNcXClcbmF4b2xvdGxfYWxsID0gcmVhZFJEUyhmaWxlID0gXFxkYXRhL2V4cHJlc3Npb24vYXhvbG90bC9wYWxsaXVtcmVzMDdfaGFybW9ueS5SRFNcXClcbmBgYFxuYGBgIn0= -->

```r
```r
axolotl_gaba = readRDS(file = \data/expression/axolotl/pallium_neuronal_GABA_res07_harmony.RDS\)
axolotl_glut = readRDS(file = \data/expression/axolotl/pallium_neuronal_Glut_res05_harmony.RDS\)
axolotl_neurons = readRDS(file = \data/expression/axolotl/pallium_neuronal_res07_harmony.RDS\)
axolotl_all = readRDS(file = \data/expression/axolotl/palliumres07_harmony.RDS\)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Load necessary eggNOG data


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuZ2FiYV9jbCA9IHBhc3RlMChcXEdBQkFfXFwsIGF4b2xvdGxfZ2FiYSRzZXVyYXRfY2x1c3RlcnMpXG5uYW1lcyhnYWJhX2NsKSA9IGNvbG5hbWVzKGF4b2xvdGxfZ2FiYSlcbmdsdXRfY2wgPSBwYXN0ZTAoXFxHbHV0X1xcLCBheG9sb3RsX2dsdXQkc2V1cmF0X2NsdXN0ZXJzKVxubmFtZXMoZ2x1dF9jbCkgPSBjb2xuYW1lcyhheG9sb3RsX2dsdXQpXG5vdGhlciA9IGF4b2xvdGxfYWxsJGNsYXNzaWZpY2F0aW9uXzJbIShjb2xuYW1lcyhheG9sb3RsX2FsbCkgJWluJSBjKG5hbWVzKGdsdXRfY2wpLCBuYW1lcyhnYWJhX2NsKSkpXVxuYWxsX2Fubm90ID0gZGF0YS5mcmFtZShcXGFsbF9hbm5vdFxcID0gYyhvdGhlciwgZ2FiYV9jbCwgZ2x1dF9jbCkpXG5heG9sb3RsX2FsbCA9IEFkZE1ldGFEYXRhKGF4b2xvdGxfYWxsLCBtZXRhZGF0YSA9IGFsbF9hbm5vdClcbmF4b2xvdGxfYWxsX2ZpbHQgPSBheG9sb3RsX2FsbFssIWlzLm5hKGF4b2xvdGxfYWxsJGFsbF9hbm5vdCkgJiAhZ3JlcGwoXFxuZXVyb25cXCwgYXhvbG90bF9hbGwkYWxsX2Fubm90KV1cbmBgYFxuYGBgIn0= -->

```r
```r
gaba_cl = paste0(\GABA_\, axolotl_gaba$seurat_clusters)
names(gaba_cl) = colnames(axolotl_gaba)
glut_cl = paste0(\Glut_\, axolotl_glut$seurat_clusters)
names(glut_cl) = colnames(axolotl_glut)
other = axolotl_all$classification_2[!(colnames(axolotl_all) %in% c(names(glut_cl), names(gaba_cl)))]
all_annot = data.frame(\all_annot\ = c(other, gaba_cl, glut_cl))
axolotl_all = AddMetaData(axolotl_all, metadata = all_annot)
axolotl_all_filt = axolotl_all[,!is.na(axolotl_all$all_annot) & !grepl(\neuron\, axolotl_all$all_annot)]

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->




# Prepare datasets
Reformat matrices to have only one2one orthologs that are shared


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxubWV0YV9saXN0ID0gbGlzdChcXGV4cF90dXJcXCA9IHR1cnRsZV9hbGxAbWV0YS5kYXRhLCBcXGV4cF9saXpcXCA9IGxpemFyZF9hbGxAbWV0YS5kYXRhLCBcbiAgICAgICAgICAgICAgICAgXFxleHBfZHJlXFwgPSBkcmVyaW9fYnJhaW5AbWV0YS5kYXRhJENsdXN0ZXIyLjUsIFxcZXhwX3plZlxcID0gZmluY2hlc19saXN0JFpGQG1ldGEuZGF0YSwgXG4gICAgICAgICAgICAgICAgIFxcZXhwX2JlZlxcID0gZmluY2hlc19saXN0JFpGQG1ldGEuZGF0YSwgXFxleHBfbW91XFwgPSBtb3VzZV9icmFpbkBtZXRhLmRhdGEsIFxuICAgICAgICAgICAgICAgICBcXGV4cF9heG9cXCA9IGF4b2xvdGxfYWxsX2ZpbHRAbWV0YS5kYXRhKVxuY2xfaWRfbGlzdCA9IGxpc3QoXFxjbHVzdGVyXFwsIFxcY2x1c3RlcnNcXCwgXFxNb2QyLjVcXCwgXFxzZXVyYXRfY2x1c3RlcnNcXCwgXFxzZXVyYXRfY2x1c3RlcnNcXCwgXFxUYXhvbm9teV9ncm91cFxcLFxuICAgICAgICAgICAgICAgICAgXFxhbGxfYW5ub3RcXClcbm5hbWVzKGNsX2lkX2xpc3QpID0gbmFtZXMobWV0YV9saXN0KVxuYGBgXG5gYGAifQ== -->

```r
```r
meta_list = list(\exp_tur\ = turtle_all@meta.data, \exp_liz\ = lizard_all@meta.data, 
                 \exp_dre\ = drerio_brain@meta.data$Cluster2.5, \exp_zef\ = finches_list$ZF@meta.data, 
                 \exp_bef\ = finches_list$ZF@meta.data, \exp_mou\ = mouse_brain@meta.data, 
                 \exp_axo\ = axolotl_all_filt@meta.data)
cl_id_list = list(\cluster\, \clusters\, \Mod2.5\, \seurat_clusters\, \seurat_clusters\, \Taxonomy_group\,
                  \all_annot\)
names(cl_id_list) = names(meta_list)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Renormalise, get mean expression and markers


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuc2V1cmF0X29uZTJvbmUgPSBsaXN0KClcbmF2Z19leHAgPSBsaXN0KClcbm1rX2xpc3QgPSBsaXN0KClcbmZvcihuIGluIG5hbWVzKGNvdW50c19saXN0KSl7XG4gIHByaW50KG4pXG4gIHNldXJhdF9vbmUyb25lW1tuXV0gPSBDcmVhdGVTZXVyYXRPYmplY3QoY291bnRzX2xpc3RbW25dXSwgbWV0YS5kYXRhID0gbWV0YV9saXN0W1tuXV0pXG4gIElkZW50cyhzZXVyYXRfb25lMm9uZVtbbl1dKSA9IHNldXJhdF9vbmUyb25lW1tuXV1AbWV0YS5kYXRhWyxjbF9pZF9saXN0W1tuXV1dXG4gIHNldXJhdF9vbmUyb25lW1tuXV0gPSBzdXBwcmVzc1dhcm5pbmdzKFNDVHJhbnNmb3JtKHNldXJhdF9vbmUyb25lW1tuXV0sIGRvLmNvcnJlY3QudW1pID0gVCwgdmVyYm9zZSA9IEYsIFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzZWVkLnVzZSA9IDEsIHZhcnMudG8ucmVncmVzcyA9IFwibkNvdW50X1JOQVwiLFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB2YXJpYWJsZS5mZWF0dXJlcy5ydi50aCA9IDEsIHJldHVybi5vbmx5LnZhci5nZW5lcyA9IEYsXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHZhcmlhYmxlLmZlYXR1cmVzLm4gPSBOVUxMKSlcbiAgYXZnX2V4cFtbbl1dID0gQXZlcmFnZUV4cHJlc3Npb24oc2V1cmF0X29uZTJvbmVbW25dXSwgYXNzYXlzID0gXCJTQ1RcIikkU0NUXG4gIG1rX2xpc3RbW25dXSA9IEZpbmRBbGxNYXJrZXJzKHNldXJhdF9vbmUyb25lW1tuXV0sIGxvZ2ZjLnRocmVzaG9sZCA9IDAuMiwgb25seS5wb3MgPSBULFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBwc2V1ZG9jb3VudC51c2UgPSAwLjEsIGFzc2F5ID0gXCJTQ1RcIilcbn1cbnNhdmUoc2V1cmF0X29uZTJvbmUsIGF2Z19leHAsIG1rX2xpc3QsIGZpbGUgPSBcInJlc3VsdHMvY3Jvc3Nfc3BlY2llc19jb3JyZWxhdGlvbi9vbmUyb25lX3NldXJhdF9hdmdfbWsuUkRhdGFcIilcbmBgYCJ9 -->

```r
seurat_one2one = list()
avg_exp = list()
mk_list = list()
for(n in names(counts_list)){
  print(n)
  seurat_one2one[[n]] = CreateSeuratObject(counts_list[[n]], meta.data = meta_list[[n]])
  Idents(seurat_one2one[[n]]) = seurat_one2one[[n]]@meta.data[,cl_id_list[[n]]]
  seurat_one2one[[n]] = suppressWarnings(SCTransform(seurat_one2one[[n]], do.correct.umi = T, verbose = F, 
                                                     seed.use = 1, vars.to.regress = "nCount_RNA",
                                                     variable.features.rv.th = 1, return.only.var.genes = F,
                                                     variable.features.n = NULL))
  avg_exp[[n]] = AverageExpression(seurat_one2one[[n]], assays = "SCT")$SCT
  mk_list[[n]] = FindAllMarkers(seurat_one2one[[n]], logfc.threshold = 0.2, only.pos = T,
                                pseudocount.use = 0.1, assay = "SCT")
}
save(seurat_one2one, avg_exp, mk_list, file = "results/cross_species_correlation/one2one_seurat_avg_mk.RData")

Correlation analysis

Load data

load("results/cross_species_correlation/one2one_seurat_avg_mk.RData")

Select genes to use for correlations, normalise data

load("results/cross_species_correlation/one2one_seurat_avg_mk.RData")

upset(fromList(top_mk_list[c(6,4,5,2,1,7,3)]), sets = names(top_mk_list)[c(6,4,5,2,1,7,3)], 
      order.by = "freq", keep.order = TRUE, number.angles = 0, nintersects = 30, point.size = 1.8)

Calculate correlations

Plot correlation matrices

sp_cor_l = list()
for(i in names(avg_exp_sub)){
  for(j in names(avg_exp_sub)){
    print(paste(i, j))
    sp_cor_l[[paste0(i, ".", j)]] = cor(avg_exp_sub[[i]], avg_exp_sub[[j]], method = "sp")
  }
}

[1] "exp_tur exp_tur"
[1] "exp_tur exp_liz"
[1] "exp_tur exp_dre"
[1] "exp_tur exp_zef"
[1] "exp_tur exp_bef"
[1] "exp_tur exp_mou"
[1] "exp_tur exp_axo"
[1] "exp_liz exp_tur"
[1] "exp_liz exp_liz"
[1] "exp_liz exp_dre"
[1] "exp_liz exp_zef"
[1] "exp_liz exp_bef"
[1] "exp_liz exp_mou"
[1] "exp_liz exp_axo"
[1] "exp_dre exp_tur"
[1] "exp_dre exp_liz"
[1] "exp_dre exp_dre"
[1] "exp_dre exp_zef"
[1] "exp_dre exp_bef"
[1] "exp_dre exp_mou"
[1] "exp_dre exp_axo"
[1] "exp_zef exp_tur"
[1] "exp_zef exp_liz"
[1] "exp_zef exp_dre"
[1] "exp_zef exp_zef"
[1] "exp_zef exp_bef"
[1] "exp_zef exp_mou"
[1] "exp_zef exp_axo"
[1] "exp_bef exp_tur"
[1] "exp_bef exp_liz"
[1] "exp_bef exp_dre"
[1] "exp_bef exp_zef"
[1] "exp_bef exp_bef"
[1] "exp_bef exp_mou"
[1] "exp_bef exp_axo"
[1] "exp_mou exp_tur"
[1] "exp_mou exp_liz"
[1] "exp_mou exp_dre"
[1] "exp_mou exp_zef"
[1] "exp_mou exp_bef"
[1] "exp_mou exp_mou"
[1] "exp_mou exp_axo"
[1] "exp_axo exp_tur"
[1] "exp_axo exp_liz"
[1] "exp_axo exp_dre"
[1] "exp_axo exp_zef"
[1] "exp_axo exp_bef"
[1] "exp_axo exp_mou"
[1] "exp_axo exp_axo"

UMAP from mean expression

cols = colorRampPalette(c(rev(RColorBrewer::brewer.pal(9, "Blues")),
                          RColorBrewer::brewer.pal(9, "Reds")))(100)
for(n in names(sp_cor_l)){
  pltmat = t(sp_cor_l[[n]])
  br = c(seq(min(pltmat), 0, length.out = 51), seq(0, max(pltmat), length.out = 51)[-1])
  pheatmap::pheatmap(pltmat, breaks = br, color = cols,
                     clustering_method = "ward.D2", main = n)
}

Network in MDS from correlations

cor_tab_l = list()
for(n in names(sp_cor_l)){
  dat = reshape2::melt(sp_cor_l[[n]])
  dat$sp1 = strsplit(strsplit(n, ".", fixed = T)[[1]][1], "_", fixed = T)[[1]][2]
  dat$sp2 = strsplit(strsplit(n, ".", fixed = T)[[1]][2], "_", fixed = T)[[1]][2]
  
  dat$Var1 = paste0(dat$sp1, "_", dat$Var1)
  dat$Var2 = paste0(dat$sp2, "_", dat$Var2)
  cor_tab_l[[n]] = dat
}
cor_tab = Reduce(rbind, cor_tab_l)
cor_sp = reshape2::dcast(cor_tab, Var1 ~ Var2, value.var = "value")
rownames(cor_sp) = cor_sp[,1]
cor_sp = cor_sp[,-1]
adj_mat = cor_sp>=0.3

# build graph, project with MDS
network_sp = graph_from_adjacency_matrix(adj_mat, weighted=T, mode="undirected", diag=F)
l_sp = igraph::layout_with_mds(network_sp)
l_sp = data.frame(l_sp)
l_sp$cl = colnames(adj_mat)
rownames(l_sp) = colnames(adj_mat)
l_sp$sp = unlist(lapply(strsplit(l_sp$cl, "_"), function(x) x[1]))
l_sp$cl = unlist(lapply(strsplit(l_sp$cl, "_"), function(x) x[2]))

df_names = l_sp[l_sp$sp=="axo",]
ggplot(l_sp, aes(x = X1, y = X2, colour = sp))+
  geom_point()+
  ggrepel::geom_text_repel(data = df_names, mapping = aes(x = X1, y = X2, label = cl), max.overlaps = 100)+
  theme_classic()

common_genes = Reduce(intersect, top_mk_list)
any_genes = Reduce(union, top_mk_list)
non_shared_genes = any_genes[!(any_genes %in% common_genes)]
non_shared_genes = non_shared_genes[non_shared_genes %in% Reduce(intersect, lapply(avg_exp, rownames))]
avg_exp_sub = lapply(avg_exp, function(x) t(apply(x[non_shared_genes,], 1, function(y) y/mean(y))))
lapply(avg_exp_sub, dim)

for(n in names(avg_exp_sub)){
  sp = strsplit(n, "_")[[1]][2]
  colnames(avg_exp_sub[[n]]) = paste0(sp, "..", colnames(avg_exp_sub[[n]]))
}
avg_exp_all = t(Reduce(cbind, avg_exp_sub))
set.seed(2954)
l = uwot::umap(avg_exp_all, metric = "cosine", ret_nn = T, n_epochs = 1000)
l_sp = data.frame(l$embedding)
rownames(l_sp) = rownames(avg_exp_all)
l_sp$cl = rownames(avg_exp_all)
l_sp$sp = unlist(lapply(strsplit(l_sp$cl, "..", fixed = T), function(x) x[1]))
l_sp$cl = unlist(lapply(strsplit(l_sp$cl, "..", fixed = T), function(x) x[2]))

df_names = l_sp[l_sp$sp=="axo",]
ggplot(l_sp, aes(x = X1, y = X2, colour = sp))+
  geom_point()+
  ggrepel::geom_text_repel(data = df_names, mapping = aes(x = X1, y = X2, label = cl), 
                           max.overlaps = 100, size = 3)+
  theme_classic()